#Clean-up and load library:
rm(list=ls())
library(mgcv)

#Load and subset data:
###SET THE WORKING DIRECTORY TO WHEREVER YOU SAVE THE DATA:###
setwd('/Volumes/MONOGAN/midterm10/variables/')
data<-read.csv('election2010_2.csv')

#Fit the estimated propensity function by modeling the treatment:
mod.treatment<-glm(treatment~repInc+demInc+propDem+unem09+forBorn+relFunds
+repInc*propDem+demInc*propDem+repInc*unem09+demInc*unem09+repInc*forBorn+demInc*forBorn+repInc*relFunds+demInc*relFunds+propDem*unem09+propDem*forBorn+propDem*relFunds+unem09*forBorn+unem09*relFunds+forBorn*relFunds, family=gaussian(link='probit'), data=data); summary(mod.treatment)

#Propensity scores are recorded in the following vector:
mod.treatment$fitted.values
#I have saved these into the variable "propensity".

#Assess balance by comparing t-ratios of models of covariates:
ideol.without<-lm(propDem~treatment, data=data);summary(ideol.without)
ideol.with<-lm(propDem~treatment+propensity, data=data);summary(ideol.with)
unem.without<-lm(unem09~treatment, data=data);summary(unem.without)
unem.with<-lm(unem09~treatment+propensity, data=data);summary(unem.with)
for.without<-lm(forBorn~treatment, data=data);summary(for.without)
for.with<-lm(forBorn~treatment+propensity, data=data);summary(for.with)
funds.without<-lm(relFunds~treatment, data=data);summary(funds.without)
funds.with<-lm(relFunds~treatment+propensity, data=data);summary(funds.with)
rInc.without<-glm(repInc~treatment, family=binomial(link=probit), data=data);summary(rInc.without)
rInc.with<-glm(repInc~treatment+propensity, family=binomial(link=probit), data=data);summary(rInc.with)
dInc.without<-glm(demInc~treatment, family=binomial(link=probit), data=data);summary(dInc.without)
dInc.with<-glm(demInc~treatment+propensity, family=binomial(link=probit), data=data);summary(dInc.with)

#Estimate the model of vote shares based on treatment effect:
outcome.model<-gam(cRepTwo~s(propensity)+s(propensity, by=treatment)+repInc+demInc+propDem+unem09+forBorn+relFunds,
           data=data, family=gaussian(link='probit'))
summary(outcome.model)

#Alternate, misspecfied model that manufactures a significant treatment effect:
bad.model<-lm(cRepTwo~treatment+repInc+demInc+unem09+forBorn+relFunds, data=data); summary(bad.model)

####CODE FOR FIGURES####

#Figure of word distribution (Figure 1 in paper):
Enforcement<-c(2.97,3.51)
Material<-c(2.41,2.28)
Security<-c(5.21,4.68)
Latinos<-c(0.22,0.17)
Ethics<-c(1.30,1.30)
mean.use<-cbind(Enforcement,Material,Security,Latinos,Ethics)
rownames(mean.use)<-c("2006","2010")
barplot(mean.use, beside=TRUE, col=c('gray20','gray80'), ylab="Mean Percentage in Category", xpd=FALSE, axis.lty=1)
legend(x=12,y=4,c(2006,2010), fill=c('gray20','gray80'))
abline(h=0, col='gray66')

#Figure of smoothed intercept term (not reported in paper):
plot(outcome.model, select=1, xlab="Propensity for Harsh Rhetoric", ylab="Smoothed Intercept", se=TRUE)

#Figure of smoothed treatment effect (Figure 2 in paper):
plot(outcome.model, select=2, xlab="Propensity for Harsh Rhetoric", ylab="Causal Effect of Immigration Rhetoric on Vote Share", se=TRUE)
